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Abstract 



o 

Advances in solid-state technology have enabled the development of silicon photomultiplier 
sensor arrays capable of sensing individual photons. Combined with high-frequency time-to- 
D digital converters (TDCs), this technology opens up the prospect of sensors capable of recording 

with high accuracy both the time and location of each detected photon. Such a capability could 
^-H lead to significant improvements in imaging accuracy, especially for applications operating with 

^-H low photon fluxes such as LiDAR and positron emission tomography. 

The demands placed on on-chip readout circuitry imposes stringent trade-off's between fill 
^ factor and spatio-temporal resolution, causing many contemporary designs to severely under- 

HH utilize the technology's full potential. Concentrating on the low photon flux setting, this paper 

t/3 leverages results from group testing and proposes an architecture for a highly efficient readout of 

, pixels using only a small number of TDCs, thereby also reducing both cost and power consump- 

tion. The design relies on a multiplexing technique based on binary interconnection matrices. 
^-H We provide optimized instances of these matrices for various sensor parameters and give explicit 

^ upper and lower bounds on the number of TDCs required to uniquely decode a given maximum 

number of simultaneous photon arrivals. 

To illustrate the strength of the proposed architecture, we note a typical digitization result 
of a 120 X 120 photodiode sensor on a 30 /xm x 30 /xm pitch with a 40 ps time resolution 
and an estimated fill factor of approximately 70%, using only 161 TDCs. The design guarantees 
registration and recovery of up to s = 4 simultaneous photon arrivals. A fast decoding algorithm 
is available, which decodes successfully with probability 1 — whenever s > 4. By contrast, 
^-H a cross-strip design requires 240 TDCs and cannot uniquely decode any simultaneous photon 

^ arrivals. In a series of realistic simulations of scintillation events in clinical positron emission 

tomography the design was able to recover the spatio-temporal location of 98.6% of all photons 
that caused pixel firings. 



1 Introduction 



Photon detection has become an essential measurement technique in science and engineering. In- 
deed, applications such as positron emission tomography (PET), single- photon emission computed 
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tomography, flow cytometry, LiDAR, fluorescence detection, confocal microscopy, and radiation de- 
tection all rely on accurate measurement of photon fluxes. Traditionally, the preferred measurement 
device has been the photomultiplicr tube (PMT), which is a high-gain, low-noise photon detector 
with a high-frequency response and a large area of collection. In particular, it behaves as an ideal 
current generator for steady photon fluxes, making it suitable for use in applications in astronomy 
and medical imaging, among others. However, they are bulky, require manual assembly steps, and 
have limited spatial resolution for PET. For these reasons, extensive research has focused on finding 
feasible solid-state alternatives that can operate at much lower voltages, are immune to magnetic 
effects, have increased efflciency, and are smaller in physical size [43]. Recent designs have raised 
significant interest in the community, and their use as a replacement for PMTs in applications such 
as PET imaging [31], high-energy physics [48], astrophysics [52], LiDAR [3], and flow cytometry 
has been recently explored. 

Silicon photomultiplier (SiPM) devices consist of two-dimensional arrays of Geiger avalanche 
photodiodes (APD) that are run above their breakdown voltage and are integrated with either 
an active or passive quenching circuit. Devices built up from these Geiger APD microcells are 
characterized by the fraction of sensing area on the device (the fill factor), and the fraction of 
incident photons that cause charged APDs to fire (the detection quantum efficiency). The product 
of these two quantities gives the combined photon detection efficiency (PDE). 

The compatibility of Geiger APDs with standard CMOS technology has enabled a number of 
different designs. For example, the digital silicon photomultiplier (DSiPM) [27] adds processing 
circuitry to count the total number of photons hitting the sensor and records the time stamp of 
the first group of photons that crosses a predetermined photon intensity threshold. The resolution 
of the resulting sensor coincides with the size of the entire sensor (usually 3 mm x 3 mm in area) 
whereas the temporal sampling is limited to a single time stamp per pulse. Nevertheless, the 
sensor does achieve a very high fill factor (80%) and is, therefore, very sensitive. Note that in 
SiPM terminology, the above tiling of microcells into larger atomic units is often referred to as a 
pixel. Throughout this paper, however, we use the term pixel to refer to individual Geiger APD 
microcells, since they represent the smallest resolvable element in our proposed design. 

An alternative sensor design, called SPAD, aims at a high temporal resolution and registers 
the time of each pixel (APD) firing [34, 40]. This is achieved by connecting each pixel to a high- 
performance time-to-digital converter (TDC), which records a time stamp in a memory buffer 
whenever a signal is detected on its input. Because of their relatively low complexity, especially 
when compared to analog-to-digital converters, TDCs allow the sensor to achieve an extremely 
high temporal resolution. However, the spatial resolution of the sensor is severely compromised by 
the large amount of support circuitry between neighboring pixels, resulting in an extremely low fill 
factor of approximately 1-5%. 

Although promising, both designs show that current implementations highly underutilize the 
full potential of these silicon devices: due to the restricted chip area a trade-off must be made 
between the spatio-temporal resolution and the fill factor. In order to improve this trade-off, we 
take advantage of the special properties present in settings with a low photon fiux. In particular, 
we claim that, by taking advantage of the temporal sparsity of the photon arrivals, we can increase 
the spatio-temporal resolution of the SiPM to within a fraction of the theoretical maximum while 
maintaining simple circuitry and a high fill factor. 

This claim is made possible by the combination of two ideas: (1) to use TDCs as the main 
readout devices, and (2) to exploit ideas from group testing to reduce the number of these devices. 
As a motivational example, consider three possible designs for an mxm pixel grid. The first design, 
illustrated in Figure 1(a) corresponds to a design with a single TDC per pixel. This design can be 
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Figure 1: Three SiPM designs with (a) one TDC per pixel, (b) one TDC per row and column of 
pixels, and (c) one TDC per bit of the binary representation of each pixel number. The first design 
is capable of detecting any number of simultaneous hits, at the cost of a large number of TDCs. 
The other two can only uniquely decode up to a single pixel firing, but require substantially fewer 
TDCs. 



seen as a trivial group test capable of detecting an arbitrary number of simultaneous firings, but 
at the cost of n? TDCs. When the photon flux is low, only a small number of pixels will typically 
fire at the same time^, thus causing most TDCs to be idle most of the time, and resulting in a very 
poor usage of resources. The second design, shown in Figure 1(b), corresponds to the widely- used 
cross-strip architecture in which rows and columns of pixels are each aggregated into a single signal. 
This design requires only 2n TDCs, but information from the TDCs can only be uniquely decoded 
if no more than one pixel fired during the same time interval. Hinting at the power of more efficient 
group-testing designs. Figure 1(c) shows a design in which pixels are numbered from 1 to and in 
which a pixel i is connected to TDC j only if the j-th bit in the binary representation number i is 
one. This design is also capable of decoding up to a single pixel firing, but requires only [log2 2nJ 
TDCs. The objective of this paper is to show that whenever the number of simultaneous firings is 
small, a significant reduction in the number of TDCs (and accompanying memory buffers) can be 
attained by carefully selected designs, which reduce both the amount of overhead circuitry and the 
generated volume of data. 

1.1 Contributions and paper organization 

The main contribution of this paper are: (1) the introduction of group testing to the design of 
interconnection networks in imaging sensors, combined with TDC readout, (2) an architecture to 
time and locate photon arrivals, (3) the construction of highly optimized group testing matrices for 
a variety of conventional m x m grid sizes, (4) a comparison with alternative constructions, and 
(5) an extensive comparison with explicitly evaluated and optimized theoretical upper and lower 
bounds on the minimum number of TDCs required to decode an imaging sensor array. 

The paper is organized as follows. We start by reviewing the concept of group testing and 
derive the currently best-known group testing designs for a range of grid sizes. We then provide an 
extensive survey of constructions used to obtain upper and lower bound on the minimum number of 
TDCs needed to guarantee recovery for d = 2, . . . , 6 simultaneous hits. There, the results show that 
the designs generated earlier are close to the theoretical lower bounds. We then test the performance 
of the design on simulated data generated using a realistic model of scintillation events arising in 

^Throughout the paper, firing at the same time is to be understood as within the same TDC samphng interval. 
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PET scanners. We conclude with a discussion on practical considerations in the implementation of 

the proposed design. 

Throughout the paper we use the following notation: log is base e unless otherwise indicated; 
[n] denotes the set {!,... ,n}; and (^^'') denotes the family of all subsets of [n] with k elements. 

2 Group testing 

Group testing was proposed by Dorfman [16] to effectively screen large numbers of blood samples for 
rare diseases (see [17] for more historical background). Instead of testing each sample individually, 
carefully chosen subsets of samples are pooled together and tested as a whole. If the test is negative, 
we immediately know that none of the samples in the pool is positive, saving a potentially large 
number of tests. Otherwise, one or more samples is positive and additional tests are needed to 
determine exactly which ones. 

Since its introduction, group testing has been used in a variety of different applications including 
quality control and product testing [50], pattern matching [11, 33], DNA library screening [4, 18, 
20, 41], multi-access communication [53], tracking of hot items in databases [13], and many others 
[17, 42]. Depending on the application, group testing can be applied in an adaptive fashion in 
which tests are designed based on the outcome of previous tests, and in a nonadaptive fashion in 
which the tests are fixed a priori. Additional variations include schemes that provide robustness 
against test errors [37, 38], or the presence of inhibitors [14, 26], which cause false positives and 
negatives, respectively. 

In our SiPM application, each pixel fires independently^ according to some Poisson process, 
and exactly fits in the probabilistic group-testing model [4]. Nevertheless, since we are interested 
in guaranteeing performance up to a certain level of simultaneous firings, we will study the appli- 
cation from a combinatorial group-testing perspective. Furthermore, measurements are necessarily 
nonadaptive as they are hardwired in the actual implementation. No error correction is needed; the 
only errors we can expect are spurious pixel firings (dark counts) or afterpulses, but these appear 
indistinguishable from real firings, and cannot be corrected for. The rate of these spurious pixel 
firings is usually much less than the signal rate, especially in our scintillation example [7]. 

2.1 Matrix representation and guEirantees 

A group test can be represented as a binary t x n incidence matrix or code A with t the length 
or number of tests and n the size or number of items, or pixels in our case. An entry Aij is set 
to one if item j is pooled in test i, that is, if pixel j is connected to TDC i, and zero otherwise. 
Given a vector x of length n with xj = 1 if item j is positive and zero otherwise, we can define 
the test vector y as y = Ax, where multiplication and addition are defined as logic AND and OR, 
respectively'^. The columns aj in A are called codewords, and for sake of convenience we allow set 
operations to be applied to these binary codewords, acting on or defining their support. As an 
illustration of the above, consider the following example: 

^This assumes that pixels have been shielded to avoid cross-talk. 

^Addition using the logic OR gives 0-1-0 = and 0-1-1 = 1-1-1 = 1, while • = • 1 = and 1 • 1 = 1, as usual. 
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In this example pixels 1, 4, and 5 are connected to the first TDC, as represented by the first row of 
A. When pixels 2 and 4 fire, we sum the corresponding columns and find that the first and second 
TDCs are activated, as indicated by y. 

For group testing to be eflFective we want to have far fewer tests than items (t ^ n). This 
inherently means that not all vectors x can be reconstructed from y, so we will be interested in 
conditions that guarantee recovery when x has up to d positive entries. 

The weakest notion for recovery is d-separability, which means that no combination of exactly 
d codewords can be expressed using any other combination of d codewords. A matrix is said to 
be d-separable if the combination of any up to d codewords is unique. This immediately gives an 
information-theoretic lower bound on the length t: 



(1) 



When d is small compared to n, this gives t k. dlog2{n/d). A stronger criteria is given by d- 
disjunctness. Given any two codewords u and v, we say that u covers v if uUv = u. Based on this, 
define A to be d-disjunct if no codeword in A is covered by the union of any d others. The concept 

of disjunctness of sets and codes has been proposed in a number of settings and such codes arc also 
known as superimposed codes [36], or d-cover free families [24]. The advantage of d-disjunct codes 
is that the positive entries in any d-sparse vector x correspond exactly to those codewords in A 
that remain after discarding all codewords that have a one in a position where y is zero. This is 
unlike general separable matrices, where a combinatorial search may be needed for the decoding. 

The central goal in group testing is to determine the smallest length t for which a matrix of 
size n can satisfy certain separability or disjunctness properties. Denote by S(n,d) and D{n,d) 
the set of all t x n matrices A that are d-separable, respectively d-disjunct. Then we can define 
Ts{n, d) the smallest length of any A G S{n, d), and likewise for Tu(n, d). As a further refinement 
in classification we define Du,{n,d) as all constant- weight d-disjunct matrices of size n, i.e., those 
matrices whose columns all have the same weight or number of nonzero entries (note that the 
weights of any two matrices within this class can differ). The overlap between any two columns 



ai and Oj of a code is defined as 



pairwise overlap 



|ai n aj|. Consider a constant- weight w code A with maximum 

H = IJi{A) := max joj H Ojj. 



It is easily seen that it takes {w/n] columns to cover another one, and therefore that 



d > 



w 



I 1^ 



(2) 



We can then define Dy,^n{n, d) as the class of constant- weight codes of size n where the right-hand 
side of (2) equals d. As for all other classes, the actual disjunctness of the codes in -0,^,^ may still 
be higher, however, there is a subtle distinction in parameterization here in that, by definition, we 
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have Dyj^n{n, d + 1) fl Dw^^{n, d) = 0, whereas Dyj{n, d + 1) C -D«,(n, d), and hkewise for D and S. 
Summarizing, we have 

Da;.i_i{n,d) C Du,{n,d) C D(n,d) C S{n,d), 
along with the associated minimum lengths T{n,d). 

2.2 Matrix construction 

In this section we discuss three methods for creating d-disjunct binary matrices. In addition, we 
construct matrices for the particular case where n = 3,600, corresponding to a 60 x 60 pixel array. 

2.2.1 Greedy approach 

The greedy approach generates d-disjunct matrices one column at a time. For the construction of a 
const ant- weight w matrix of length t the algorithms starts with J-" = 0, and proceeds by repeatedly 
drawing, uniformly at random, an element F from as a candidate column. Whenever the 
distance to all sets already in T exceeds w/dwe add F to T and continue the algorithm with the 
updated family. This process continues until either | = n or a given amount of time has passed. 

We applied the greedy algorithm to construct d-disjunct matrices of size at least 3,600. For 
each value of d and length t, the algorithm was allowed to run for 12 hours. Row (1) of Table 3 
gives the minimum t for which the algorithm found a solution. Instances with fewer rows are very 
well possible; this strongly depend both on the choice of the initial few columns and the amount 
of time available. 

A further reduction in the code length could be achieved by checking disjunctness in a groupwise, 
rather than a pairwise setting. However, the number of possible d-subsets of columns to consider 
at each iteration grows exponentially in d, thereby rendering this approach intractable for even 
moderate values of d. 

2.2.2 Chinese remainder sieve 

Eppstein et al. [22] recently proposed the 'Chinese remainder sieve', a number-theoretic method 
for the deterministic construction of d-disjunct matrices. In order to generate a d-disjunct matrix 
with at least n columns, first choose a sequence {pl^,P2^, ■ ■ ■ ,Pk'°} of powers of distinct primes such 
that HiLiPi' > n'^- At X n matrix A with t = Yli=iPT then created by vertical concatenation 
of X n matrices as follows: 



As an example, consider the construction of a 1-disjunct 5x6 matrix generated based on pi = 2 
and p2 = 3. Applying the above definition we obtain: 
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d t n prime powers 

^ 82 4077 {2'^3^5,7, 11,13, 17,23} 

3 155 3855 {22,3,5,7,11,13,17,19,23,29,37} 

4 237 3631 {2^,3,5,7,11,13,17,19,23,29,31,37,43} 

5 333 4023 {23,32,5,7,11,13,17,19,23,29,31,37,41,43,53} 

6 445 4077 {23,32,5,7,11,13,17,19,23,29,31,37,41,43,53,61} 

Table 1: Prime powers used for the construction of d-disjunct t x n. 

The above construction requires t = 0{d'^log^ n/{logd + log log n)) rows for a rf-disjunct matrix 
of size n [22]. This is only slightly worse than the best-known 0(d'^ log n) construction, which we 
discuss in more detail in Section 3.1. Using a near-exhaustive search over prime-power combinations, 
we construct d-disjunct matrices of size at least 3,600. The parameters giving designs with the 
smallest number of rows are given in Table 1. 

2.2.3 Designs based on error-correcting codes 

Excellent practical superimposed codes can be constructed based on error-correcting codes. Here, 
we discuss a number of techniques and constructions we used to generate the best superimposed 
codes known to us for a variety of pixel array sizes (i.e, for square m x m arrays with m G 
{10,20,30,40,60,120}). For comparison with other constructions, the results are summarized in 
Table 2, and in row (n) of Table 3 for n = 3,600. 

Binary codes. The most straightforward way of obtaining d-disjunct superimposed codes is by 
simply taking constant- weight binary error-correction codes obeying (2), with overlap /x as given 
below. The on-line repository [6] lists the best known lower bound on maximum size A(n, d, w) for 
constant weight w codes of length n and Hamming distance d (note the different use of n and d in 
this context). Given a code ST(n,d,w) from this Standard Table, the overlap satisfies n < w — d/2. 
Some codes are given explicitly, whereas others require some more work to instantiate. We discuss 
two of the most common constructions that are used to instantiate all but one of the remaining codes 
we use. The first construction consists of a short list of seed codewords Vi, along with a number 
of cyclic permutations. These permutations give rise to a permutation group V and the words in 
the final code are those in the union of orbits of the seed vectors: Ui{w | w = P{vi),P G V}. The 
second construction shortens existing codes in one of two ways. We illustrate them based on the 
generation of the 2-disjunct code St(21,8,7) of size 100 from st(24,8,8). The first type of shortening 
consists of identifying a row i with the most ones and then selecting all columns incident to this row 
and deleting row i. This both reduces the weight and the length of the code by one, but preserves 
the distance, and in this case gives St(23,8,7) (note that shortenings do not in general lead to new 
optimal or best known codes). The second type of shortening identifies a row i with the most 
zero entries and creates a new code by keeping only the columns with a zero in the i-th row, and 
then removing the row. This construction does not affect the weight of the matrix, only the size 
and length. Repeating this twice from St(23,8,7) yields the desired St(21,8,7) code. Note that for 
constant- weight codes, the minimum number of ones in any given row is bounded above by [wn/t\ . 
This expression can be used to give a theoretical minimum on the number of rows by which we can 
shorten a code. In practice it may be possible to remove a substantially larger number of rows. 
Below we will frequently use shortening to obtain codes with smaller length. In these cases we only 
use the second type of shortening, based on zero entries. 
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q-ary error-correction codes. The Standard Table only lists codes with relatively small lengths 

and consequently limited sizes. In order to construct larger or heavier codes we apply a construction 
based on (maximal distance separable (MDS)) q-arj error-correcting codes, as proposed by Kautz 
and Singleton [36]. Let {n,k,d)q denote a linear q-aiy code of length n, size q'', and Hamming 
distance d. Each codeword in these codes consists of n elements taken from GF(q') and differs in 
at least d locations from all other codewords. A binary superimposed code can be obtained from 
these g-ary codes by replacing each element with a corresponding column of a qx q identity matrix 
Iq. That is, we map each element in GF((;) to a unique column of Iq. For example, we map value 
k to the {k + l)-st column of Iq as follows: 

q-avy 2 12 110 

1 1 
binary i 1 1 1 
2 1 1 

The overlap between any two codewords of the resulting concatenated code is bounded by the 
length of the code n minus the distance d. Meanwhile, the weight is exactly the length. The 
disjunctness of the resulting code is therefore at least \{n — l)/(n — d)~\. As an aside, note that 
this construction requires an explicit set of codewords. This can be contrasted with g-ary error- 
correction codes, for which fast encoding-decoding algorithms often exist, and which do not require 
such an explicit representation. 

Existence of a large class of (n, k, d)q codes, including the well-known Reed-Solomon codes is 
shown by Singleton [36, 49]. In MacWilliams and Sloane [39, Ch. 11, Thm. 9] cyclic MDS codes 
with n = q+1 and d = q—k+2 are shown to exist whenever l<k<q+l and g is a prime power. As 
a result, it can be concluded that, when expressed in group-testing notation, this concatenated code 
construction requires a length t = C'(min[n, A;^ log^ n]) [42], compared to the best-known bound of 
0{k^\ogn). Despite the slightly weaker bound, we shall see below that for small instances, the 
resulting codes are far superior to the random constructions used to yield the improved bound. 

As an example, wc used a concatenation of (10, 4, 7)ii with /n, denoted (10, 4, 7);^^ to construct 
the 3-disjunct matrix of size n = 14,400 shown in Table 2. Most of the other codes obtained using 
this construction have an additional superscript s{k) to indicate the application of k shortening 
steps. The 4-disjunct code of length n = 1,600 also has a superscript x to indicate extension of 
the code. In this case, the four shortening steps resulted in a 113 x 1,596 code, falling just short 
of the desired size of 1,600. The very structured nature of these concatenated codes means that 
many constant-weight vectors are not included, even if they may be feasible. We can therefore try 
to apply greedy search techniques to augment the code. For this particular case it was found to be 
relatively easy to add several columns, thus resulting in a code of desired size. 

All of the g-ary codes appearing in Table 2 are Reed-Solomon codes, except for (10, 4, 7)9, which 
is a constacyclic linear code whose construction is given by [30]. 

General concatenated codes. As mentioned by Kautz and Singleton [36], it is possible to 

replace the trivial identity code by arbitrary d-disjunct binary matrices in forming concatenated 
codes, provided that the size of the matrix is at least q. This construction is extensively used 
by D'yachkov et al. [19] to form previously unknown instances of superimposed codes. We also 
investigated this approach and found that the concatenation of the (7, 4, 4)ii code with st(9,4,3) 
yielded our smallest 2-disjunct code of size 14,400. Likewise a 2-disjunct 51 x 4,489 matrix can 
be obtained by concatenating (3,2,2)67 with st(17,6,5). The code given in Table 2 lists instead 
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n = 100 n = 400 n = 900 


d 


t construction 


t construction 


t construction 


2 
3 

4 
5 
6 


21 STf21 8 7) 
36 St(36,10,7) 
48 St(48,8,5) 
60 st(60,10,6) 
75 (7,2,6){^i''^^^ 


31 STf31 8 7) 
51 St(51,10,7) 
64 0(65,9,3)^(1) 
107 (ll,3,9){f ^^^^ 
144 (13,3,11){|5'^^^^^ 


38 STf38 8 71 
73 (7,3,5);i'^^') 
95 (9,3,7){l''^^^ 
117 (ll,3,9){^i''^^^ 
156 (13,3,11){|'^^^^^ 


n = 1,600 n = 3,600 n = 14,400 


d 


t construction 


t construction 


t construction 


2 
3 
4 
5 
6 


44 st(44,8,7) 
78 (10,4,7)9"*^^^^ 
113 (9,3,7){|''^^^''' 
140 (11,3,9){|''^^^ 
166 (13,3,11){|'^^^^ 


51 St(51,8,7) 
85 (10,4,7)J^'''"'^ 
142 (9,3,7){f 
174 (ll,3,9){f('^ 
206 (13,3,11){^/^^^ 


63 (7,4,4)^,^'(^'^'^) 
110 (10,4,7){l 
161 (13,4,10){|''^®^ 
233 (16,4,13){^6''^^^^ 
323 (13,3,11)J|'*^^^ 



Table 2: Overview of best-known constructions for d-disjunct matrices of size n. Codes from the 
Standard Table [6] and covering are indicated by St(-) and c(-); all other entries are {n,k,d)q 
codes. Superscripts s{k) indicates a shortening of k steps with pivoting on zero entries; x indicates 
a greedy extension, Iq indicates concatenation with a q x q identity matrix. The St(-) superscript 
indicates concatenation with a binary code from the Standard Table. 

St(51,8,7), which yields a slightly smaller code with weight 7 instead of 15. Note that concatenation 
with non-trivial disjunct matrices can result in codes in Dy,{d,n) \ Dy]^f^{d,n) or even D{d,n) \ 
Dyj{d,n). 

Other constructions. We would like to mention two other construction techniques that can 

be used to generate superimposed codes. The first technique is based on certain designs such as 
t-designs and Steiner systems [12], which were used for example by Balding et al. [4] to construct a 
4-disjunct 65 x 520 binary code. The matrix corresponding to this particular system can be found 
on-line in the La Jolla Covering Repository [29] as C(65, 9, 3). In Table 2 we shorten this matrix 
by a single row to obtain a code for n = 400. Although this construction is significantly shorter 
than the best q-ary based design we could find (a 76 x 408 shortening of (9,3,7)9), it is typically 
much harder to obtain explicit instantiations of designs, or even to show their existence, at least 
compared to g-ary codes. The second technique based on Latin squares appears in [36], and was 
used by D'yachkov et al. [19] to create 2-disjunct 51 x 4,624 and 63 x 14,400 matrices. 

3 Theoretical bounds 

We now summarize an extensive literature concerning theoretical bounds on the minimum number 
of pools required for a group-testing design of size n to be d-disjunct. These results include both 
upper and lower bounds on this minimum, and can be used to get an idea about the quality of the 
matrices we constructed in the previous section. 
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Figure 2: Comparison between the number of TDCs used for the best-known codes in Table 2 and 
a standard cross-strip design, for different numbers of pixels n. The cross-strip design connects 
each row and column to a single TDC, giving a d = 1 disjunct encoding. The number of TDCs in 
the cross-strip design grows as 0{^/n), much faster than the 0{d? logn) in group-testing matrices. 

3.1 Asymptotic bounds 

Bounds on the growth rate of Tu(n, d) have been discovered and rediscovered in different contexts 
in information theory, combinatorics, and group testing [47]. In the context of superimposed codes, 
D'yachkov and Rykov [21] obtained the following bounds: 

n{(f\ogn/\ogd) < TD{n,d) < ©(d^logn). 

Ruszinko [47] and Fiiredi [28] give an interesting account on the lower bound and provide simpler 
proofs. The lower bound was also obtained for sufficiently large d with d < n^^'^ by Chaudhuri and 
Radhakrishnan [8] in the context of d-cover free systems for system complexity, which was extended 
to the general case by Clementi et al. [10]. 

For the upper bound, it follows from the analysis of a greedy-type construction given by Hwang 
and Sos [32], that for t > 16(i^ we have 

To^A^.d) < 16a(2 . log3(2) • (log2(n) - 1) < lldHog^{n). (3) 

An efficient polynomial-time algorithm for generating similar constant-weight d-disjunct t x n ma- 
trices with t = 0{d^\ogn) was given by Porat and Rothschild [42]. 

3.2 Non-asymptotic results 

When working with particular values of d and n, constants in theoretic bounds become crucial. 
Most theoretical results, however, are concerned with growth rates and even if explicit constants 
are given, such as those in (3), they may be too loose to be of practical use. 

In this section we are interested in the code length t required to guarantee the unique recovery 
of up to d = 6 simultaneous firings on a 60 x 60 pixel array, that is n =3,600. We do this by 
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studying the various models used to obtain upper and lower bounds on T(n, d) and numerically 

evaluating the original expressions, instead of bounding them. By optimizing over parameters such 
as weight or matrix size, we obtain the best numerical values for the bounds, as permitted by the 
different models. These will be used as a reference point to evaluate the quality of the codes given 



All bounds are evaluated using Sage [51] and summarized in Table 3 for d = 2, . . . , 6. We omit 
results for d = 1, which only requires the columns of the matrix to be unique. 

3.2.1 Upper bounds on r(n, cZ) 

The construction of a d-disjunct binary t x n matrix, or providing an algorithm for doing so, 
immediately yields t as an upper bound for one of the T(n, d) values. In the literature there are 
three main techniques to obtain d-disjunct matrices of various sizes; we discuss each of these in the 
following paragraphs. 

Sequential picking. The construction given by Hwang and Sos [32] for obtaining an upper bound 

on Tf)^ ^ (n, d) works as follows. Starting with a family of initially all good sets Go = (^') of a fixed 
weight w we pick a random element g ^ Gq. After we make this choice we determine the family Bi 
of sets in Go whose overlap with g is w/d or more (this includes g itself). We then remove all the 
bad items Bi from Go to get Gi = Go \ -Bi . We then repeat the procedure, picking at each point 
a g from Gi, determining Bi^i and forming Gj+i, until G„ is empty. We can then form a t x n 
matrix A with the support of each column corresponding to one of the selected g. This matrix will 
be d-disjunct since, by construction, it has constant column weight w and satisfies (2) with overlap 
fi < w/d. For the number of rows in A, notice that the size of each Bi satisfies 



This can be turned into the upper bound on Tp^^ (n, d) shown in (3) by following the parameter 
choice and analysis given in [32]. Essentially the same argument was used earlier by Erdos et al. 
to prove an alternative bound given in [24, Proposition 2.1]. Entries (a~c) in Table 3 are obtained 
using respectively (3), the bound given in [24] with optimal choice of w, and the smallest t that 
satisfies (4) for some w <t (see also [15, p. 232]). 

Random ensemble. An alternative to picking compatible columns one at a time, is to draw 
an entire set of columns from a suitable distribution and check if the resulting matrix is indeed 
d-disjunct. An upper bound on the required number of rows is obtained by finding the smallest t 
such that the probability that the matrix does not satisfy d-disjunctness is strictly bounded above 
by one. 

Consider a t x n matrix A with entries drawn i.i.d. from the Bernoulli distribution; each entry 
takes on the value 1 with probability P and the value with probability 1 — p. For any column in 
A, the probability that it is covered by d < n of the other columns is given by 



in Table 2. 




Because the size of the initial family Go was (^) we have that 




(4) 



p= (1-/5(1-/3)'^)*, 



(5) 
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since each entry in the column is not covered by the union of the other d if and only if its value is 
1, and the corresponding entries in the other d columns are all 0. Since we want to minimize the 
chance of overlap, we find the minimum with respect to /3, giving /3 = 1). By taking a union 

bound over all possible sets, we can bound the probability that at least one column in A is covered 
by some d others as 

P{A not d— disjunct) < ^ 

This setup forms the basis for the analysis given by D'yachkov and Rykov [21]. In Tabic 3, row 
(d), we show the smallest t for which this probability is less than one; for the optimal /3, this is 
given by 

'=['-'°K"("rf'))/'°K'"(dTT)«, 

To construct matrices with constant- weight w, wc can uniformly draw columns from S = (^). 
The probability that a fixed column is covered by d others can be evaluated as follows. Let Rd{t, w, i) 
denote the number of ways we can cover a fixed set of w entries in a t-length vector with d columns, 
given that i entries have already been covered. With only a single column left, we first cover the 
remaining w — i entries, leaving i entries to be freely selected in the remaining t — {w — i) positions. 
This gives 

jy (. .X ft-w + i^ 
Ki{t, w, I) = 

For d > 1 we can cover up to — i entries and cover the others using the remaining d—1 columns. 
Therefore, 

fw — i\ft — w+j\ 
Rd{t,w,i) = . jRd-i{t,w,t + J). 

The probability that a column is covered by d others is then given by 

p = Rd{t,w,0)/(^^^ , (7) 

and we can again use the union bound given by (6). In order to determine the smallest possible t 
we start at t = 1 and double its value until the right-hand side of (6) is less than one and then use 
binary search to find the desired t. This is repeated for all suitable column weights w and row (e) 
of Table 3 lists the smallest value of t obtained using this procedure. 

The use of the union bound in (6) ignores the fact that many of the columns and d-sets are 
independent. In order to obtain sharper bounds for a similar problem, we follow Yeh [54] and apply 
Lovasz Local Lemma [25], stated in [2, Corollary 5.1.2]. 

Lemma 3.1 (The Local Lemma; Symmetric case) Let Ei, E2,. . . ,En be events in an arbi- 
trary probability space. Suppose that each event Ei is independent of all other events Ej except at 
most fi of them, and that P{Ei) < p for all 1 < i < n. If 

ep(M + l)<l (8) 

then P(nf^i^i) > 0. 
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Working again with n columns drawn uniformly and independently from (f*'), let Ei^j denote 
the event that column i is covered by the set J of d other columns. Any event Ei^^j-^ is independent 
from Ej,^,,].-, whenever (Ji U {ii}) n (J2 U {^2}) = 0, and the number of events /x that violate this 
condition is given by 

n \ fn — d — 1^ 
^d+\) ~\d+\ 

All wc then need to do is find the smallest t for which, with an optimal choice of condition (8) 
is satisfied using p as given in (7). That is, find a t such that 

e-Rd{t,w,())(^^ {^^ + l)<l. 

The results obtained by this method are given in row (f) in Table 3. We omit similar results 
obtained using a Bernoulli model, or the g-ary construction considered by Yeh [54] , since they give 
larger values of t. 

Random selection with post processing. A third approach to generating d-disjunct matrices 
is to start with a t x m, matrix with m > n, drawn from a certain random ensemble, and randomly 
mark for deletion one of each pair of columns whose overlap is too large. Whenever the expected 
number of columns marked for deletion is strictly less than n— m+1 we can be sure that a d-disjunct 
t X n matrix exists and can be generated using this approach. (Note that some of the approaches 
described in the previous paragraph are a special case of this with m = n.) A variation on this 
approach is to pick a column and mark it for deletion if there exist any other d columns that cover 
it. Starting with this latter approach for the Bernoulli model with /3 = l/((i+l), a bound on the 
probability of a column being marked is found by applying a union bound to (5): 

P(oovered)<(''";l)(l ^^;;,,,;* 

The expected number of columns marked is then bounded by 

- 1\ / d'^ ^ * 

E(#cols) = m ■ P(covered) < m 1 



d J \ {d + iy+^ J ' 

and the smallest t for which the left-hand side of this inequality is strictly less than m — n + 1 is 
given in row (g) of Table 3. A similar derivation based on constraints on the pairwise overlap of 
columns can easily be seen to yield 

nix /t\ \ - w\ t — w 



E(#columns marked) < ■ ^ ^ 



i=\w/d\ 



W — I 



which leads to the results shown in row (h). The best results for T£,{n, d) based on bounds of the 
probabilistic method in row (i) were obtained for groupwise covering of constant weight vectors, 
a model which was studied earlier by Riccio and Colbourn [45]. For completeness we list in row 
(j) the results based on a random g-ary construction similar to the one we discussed earlier. A 
final method we would like to mention chooses m = n, but instead of removing marked columns, it 
augments the matrix with the identity matrix below all marked columns, thus fixing a (potential) 
violation of d-disjunctness by increasing the number of rows. The results of this method by Yeh [54] 
are presented in row (k). 
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Bound \ 


d 


2 


3 


4 


5 


6 


Reference and comments 


Upper bounds (theoretical) 


(a) 


Tn 


n. dj 


< 


436 


982 


1746 


2729 


3929 


Adaptation of [32] given in (3) 


(b) 


Tn 


[n, d) 


< 


94 


237 


443 


711 


1043 


[24, Proposition 2.1] with optimal w 


fc) 


Tn 


n, d) 


< 


90 


222 


412 


660 


966 


Equation (4) with optimal w 


(d) 


Tn 


n, d) 


< 


149 


278 


442 


640 


870 


Bernoulli, [21, Theorem 5] 


(e) 


Tn 


[n, d) 


< 


96 


190 


312 


459 


631 


Constant weight, groupwise cover 


(f) 




[n, d) 


< 


76 


163 


279 


422 


590 


Constant weight, Lovasz local lemma 


(g) 


Td{ 


n, w) 


< 


110 


225 


376 


561 


779 


Bernoulli, groupwise cover 


(h) 




[n, d) 


< 


99 


249 


465 


746 


1094 


Constant weight, pairwise overlap 


(i) 




[n, d) 


< 


71 


154 


265 


402 


565 


Constant weight, groupwise cover; see [45] 


(j) 


TDy,,.f, 


[n, d) 


< 


130 


300 


522 


828 


1178 


Adaptation of [9, Algorithm 1] 


(k) 


Td 


[n, d) 


< 


98 


205 


348 


524 


734 


[54, Theorem 2.6] 


Upper bounds (instances) 


(1) 




[n, d) 


< 


58 


132 


224 


345 


484 


Greedy search 


(m) 


Toy, 


[n, d) 


< 


82 


155 


237 


333 


445 


[22, Section 2]; see Table 1 


(n) 




[n. d) 


< 


51 


85 


142 


174 


206 


Error-correction code based 
















Lower 




.s 


(o) 


Toy,^^ 


[n, d) 


> 


35 


56 


75 


96 


115 


[35, Theorem 2] 


(P) 


TOyj^f, 


[n, d) 


> 


35 


55 


74 


93 


113 


[36, Equation 5] and [21, Equation 18] 


(q) 




[n, d) 


> 


35 










[23, Theorem 1] 


(r) 


TOy,,^ 


[n, d) 


> 


32 


48 


62 


75 


88 


[24, Proposition 2.1] 


(s) 


TOy, 


[n, d) 


> 


29 


40 


48 


55 


61 


[47, Lemma 3.2] 


(t) 


Ts 


[n, d) 


> 


23 


33 


43 


53 


62 


Information theoretical lower bound 



Table 3: Upper and lower bounds on the minimum number of rows needed for a d-disjunct or d- 
separable matrix with n = 3600 rows. Entries shown in boldface denote the best theoretical upper 
and lower bounds we could find for each of the given classes. 



Remarks on the upper bounds. Even though different models may lead to the same growth 
rates in T{n,d), there is a substantial difference in resulting values, when numerically evaluated. 
The groupwise models easily outperform bounds based on pairwise comparison overlap. In addition 
there is a large gap between matrices with constant-weight columns, and those generated with i.i.d. 
Bernoulli entries, the former giving much more favorable results. The second approach given above, 
based on drawing fixed t x n matrices is clearly outperformed by the third approach in which a 
t X m matrix is screened and reduced in size by removing columns that violate disjunctness or 
maximum overlap conditions, even if the Lovasz Local Lemma is used to sharpen the bounds. 

3.3 Lower bounds on T(n.d) 

Most of the lower bounds listed in Table 3 are derived using the concept of a private set. Let 
be a family of sets (the columns of our matrix A). Then T is a private set of F e T, ii T C F, 
and T is not included in any of the other sets in J^. For G Dw^^{n,d) with length t, constant 
column weight w, and maximum overlap /i, it is easily seen that, by definition, each (^u + l)-subset 
of F G J-" is private. The number of such private sets \F\ ■ cannot exceed the total number of 
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these sets, ( ^^), and it therefore follows [21, 36] that 




p + l 



w 



) 



Based on this we can find the smallest t for which there are values w and n obeying {w — 1)//U > d, 
such that the right-hand side of the above inequality exceeds the desired matrix size. This gives a 
lower bound on the required code length and the resulting values are listed in row (p) of Table 3. 
A slightly better bound on the same class of matrices follows from an elegant argument due to 
Johnson [35], and is given in row (o). 

Ruszinko [47] studies d-disjunctness in constant-weight matrices without considering maximum 
overlap and provides the following argument. Let T G Dyj{n,d) be a family of u;-subsets in [t], 
such that no d sets in cover any other set. Moreover, assume that w = kd for some integer k. 
Any F ^ T can be partitioned into d sets of length k, and it follows from the d-disjunctness of J- 
that at least one of these d subsets is private to F (otherwise it would be possible to cover F with d 
other columns). The key observation then, is that Baranyai's Theorem [5] guarantees the existence 
of s = {^i^/d different partitions of F such that no subset in the partitions is repeated. Each 
of these s partitions contains a private set, so the total number of private sets is at least \F\ ■ s. 
Rewriting this gives 



With proper rounding this can be extended to general weights giving the results shown in row 
(s). The results in rows (q) and (r) are also obtained based on private sets, but we will omit the 
exact arguments used here. Finally, an evaluation of the information-theoretic bound (1) is given 
in row (t). 

4 Numerical Experiments 

To illustrate the strength of the proposed design, we analyze the performance of some of the group- 
testing designs given in Table 2, when used to decode scintillation events in a PET setting. 

4.1 Simulation parameters 

In PET, scintillation crystals are used to convert 511 kilo-electron volt annihilation photons into 
bursts of low-energy (visible light) photons. We use the standard software package Geant4 [1, 46] to 
simulate this process for some 1,000 scintillation events for a 3 x 3 x 10 mm^ cerium-doped lutetium 
oxyorthosilicate (LSO) crystal (7.4 g/cm^, n = 1.82, absorption length = 50 m, scintillation yield = 
26,000 photons/MeV, fast time constant = 40 ns, yield ratio = 1, resolution scale = 4.41) coupled 
to a 3 X 3 X 0.75 mm silicon sensor by 50;U of optical grease (n = 1.5, absorption length = 50m). 
For each event the simulation yields the location and arrival time of the low-energy photons with 
respect to the start of the event. 

The silicon sensor is assumed to have a 70% fill factor and a quantum efficiency of 50% for blue 
photons. Once a photon is detected a pixel will be unable to detect a new photon for a known dead 
time, simulated to lie between 10 and 80 ns. We assume that each pixel takes up the same fraction 
of detector area and assign photons uniformly at random to a pixel. This increases the pixel firing 
rate by avoiding hitting dead pixels in high flux regions, thereby making the decoding process more 
challenging. We do not model dark counts and cross-talk between pixels, since they would account 
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to less than 1% of the total number of detected photons and therefore do not significantly affect 
the result. 

In terms of pixel firings, we ignore the jitter in the time between an incident photon and the 
actual firing of a pixel. This assumption greatly simplifies the simulation, and is not expected to 
have any significant influence on the results. The signal delays over the interconnection network 
between pixels and TDCs are assumed to be uniform, i.e., the travel time of the signal from each 
pixel to any connected TDC is assumed to be the identical, and can therefore be set to zero without 
loss of generality. TDCs are further assumed to be ideal in the sense that they have no down time; 
a time stamp is recorded whenever an event occurred during the sampling interval. 

4.2 Decoding 

Pixel firings give rise to specific patterns of TDC recordings through the group-testing design 
embedded in the interconnection network. These patterns can be decoded by considering the 
time stamps recorded in the memory buffer associated with each TDC and forming a binary test 
vector for each time interval. Given a test vector y, the decoding process starts by identifying 
all codewords aj that are covered by y. When the group testing matrix A is a d-disjunct matrix, 
it immediately follows that whenever there are no more than d simultaneous firings, only the 
codewords corresponding to those pixels will be selected, and the decoding is successful. Whenever 
there are s > d pixels that fired simultaneously, many more than s columns may be covered by 
y. If none of these columns can be omitted to form y then those columns coincide with pixels 
that fired and decoding is again successful. Otherwise, the decoding is ambiguous and considered 
unsuccessful. When decoding is successful we recover the pixels that fired, otherwise they are 
missed. 

4.3 Results 

In our simulations we considered sensors with both 60 x 60 and 120 X 120 pixel arrays, and a variety 
of different pixel dead times and TDC interval lengths. Due to space limitations we can only show 
a select number of tables that are representative of the results, and we will describe any significant 
differences in the text. 

Table 4 shows the simulation results obtained using 60 x 60 and 120 x 120 arrays with inter- 
connection networks based on the superimposed codes from Table 2, with disjunctness ranging 
from 2 to 6. The pixel dead time and TDC interval length are chosen to be 20 ns and 40 ps, 
respectively. These times are somewhat pessimistic; the typical dead time is longer, while the TDC 
interval could be shorter. Both choices cause an increase in the number of simultaneous photons 
per TDC sampling window, thereby making recovery more challenging. In particular, note that 
the maximum number of simultaneous firings reaches 14, far above the guaranteed recovery level. 
Nevertheless, it can be seen that even with a 4-disjunct interconnection network, more than 99% 
of all pixel firings are successfully recovered. For the 120 x 120 array, this mimbcr is slightly lower 
at 98.5% since fewer photons hit dead pixels, thus causing an increase in the number of firings and, 
consequently, an increase in the number of simultaneous hits. 

A more complete picture on the relationship between the number of pixel firings missed and 
the pixel dead time and TDC interval length is given in Tabic 5. Shown are the percentage of pixel 
firings missed for various pixel dead times and TDC interval lengths. As expected, smaller TDC 
intervals reduce the number of simultaneous hits and therefore lead to more uniquely decodable 
events. 
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(a) 60 X 60 pixel array 





Disjunctness (#TDCs) 
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(b) 120 X 120 pixel array 



Table 4: Decoding statistics for two pixel arrays using a pixel dead time of 20ns and TDC interval 
length of 40ps. Results shown correspond to the decoding of 1,000 sequential scintillation events. 
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Dead time TDC #Simult. 
(^firings) interval 


Pixel firings missed (%) 
d=2 d=3 d=A d=5 d=6 


10ns 5ps 6 
(3.4m) lOps 7 
20ps 11 
40ps 14 


1.32 0.02 0.00 0.00 0.00 
4.69 0.25 0.00 0.00 0.00 

14.48 1.87 0.04 0.00 0.00 
35.39 10.52 0.87 0.13 0.02 


20ns 5ps 6 
(3.1m) lOps 7 
20ps 11 
40ps 14 


1.21 0.02 0.00 0.00 0.00 
4.31 0.23 0.00 0.00 0.00 
13.30 1.74 0.04 0.00 0.00 
32.57 9.63 0.83 0.13 0.02 


40ns 5ps 6 
(2.8m) lOps 7 
20ps 11 
40ps 14 


1.22 0.03 0.00 0.00 0.00 
4.33 0.25 0.00 0.00 0.00 
13.21 1.83 0.05 0.00 0.00 
31.60 9.95 0.90 0.14 0.02 


80ns 5ps 6 
(2.5m) lOps 7 
20ps 11 
40ps 14 


1.34 0.03 0.00 0.00 0.00 
4.73 0.28 0.00 0.00 0.00 
14.33 2.03 0.05 0.00 0.00 
33.83 11.04 1.00 0.16 0.03 



(a) 60 X 60 pixel grid 
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(b) 120 X 120 pixel grid 



Table 5: Percentage of pixel firings missed for various combinations of pixel dead time, TDC interval 
length, and disjunctness. Results shown correspond to the decoding of 1,000 sequential scintillation 
events. The first three columns show the dead time (and total number of pixel firings), TDC time 
interval, and maximum number of simultaneous hits, respectively. 
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Figure 3: Percentage of successfully decoded vectors as a function of sparsity level for arrays of 
size 60 X 60 (solid) and 120 x 120 (dashed). For each sparsity level 10,000 random vectors were 
generated and decoded. 

As seen, even when the number of simultaneous pixel firings exceeds the disjunctness of the 
group-testing design, it often remains possible to uniquely decode the resulting codeword. To get 
more accurate statistics, we studied the decoding properties of randomly generated sparse vectors 
(corresponding to random pixel firing patterns) with the number of nonzero entries ranging from 
1 to 20. For each sparsity level we decoded 10,000 vectors and summarize the success rates in 
Figure 3. The plots show that recovery breaks down only gradually once the sparsity exceeds the 
disjunctness level of the matrix. 

5 Discussion 

For the adaptation of the group-testing based design in practical designs, a number of additional 
issues needs to be addressed. We discuss some of the major ones below. 

Decoding. Since the TDCs store time stamps in a memory buffer, we can perform the decoding 
process off-line. This is analogous to the operation of a digital storage oscilloscope. Aside from 
the obvious benefit of reducing the amount of circuitry, it has the advantage that we have a global 
view of the entire data, not just the current event, or a short history of previous events. This 
data can be used to extract more information from unsuccessfully decoded events by excluding 
columns corresponding to pixels that fired shortly before or after the current event, or by using 
more sophisticated, but perhaps combinatorial decoding techniques. 

Asynchronicity of the system. One challenge in the implementation of the proposed design is 
to make sure that the TDCs connected to a single pixel all record the firing of that pixel in the same 
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TDC interval. In other words, it is important that codewords are recorded in a single interval and 

not spread out over two consecutive intervals. Such a shift can occur due to differences of travel 
times of signals between the pixel and the TDCs as a result of differences in wiring lengths. The 
appearance of partial codewords in test vectors complicates decoding and may cause pixel firings 
to be masked by others and be missed completely with no means of detection. 

There are two main approaches of dealing with this problem. The first approach is to synchronize 
the logical signal generated by the pixel using a global clock across the chip to ensure uniform 
arrival of signals to the TDC digitizers. Such an approach has been successfully implemented in 
clock design trees of CPU chips that maintain clock skews less than the TDC sampling frequency 
[44]. The second approach is to consider bursts of consecutive events and decode the union of the 
corresponding test vectors as a whole. This also leads to an increase in the expected number of 
simultaneous firings as well as a slight decrease in temporal resolution. This approach is possible 
only when the photon flux is sufficiently low; when many consecutive TDC windows contain events, 
their union simply contains too many pixel firings to be successfully decoded. 
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